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Abstract 

In  this  paper,  the  new  higher  order  wall  boundary  conditions  are  proposed  for  solving  the 
incompressible  flows.  The  square  driven  cavity  flows  are  simulated  by  using  the  variable  order 
method  of  lines  with  the  present  wall  boundary  conditions.  The  variable  order  method  of  lines  is 
constructed  by  the  spatial  discretization,  i.e.,  the  variable  order  proper  convective  scheme  for 
convective  terms  and  the  modified  differential  quadrature  method  for  diffusive  terms,  and  time 
integration.  The  2nd,  4th,  6th,  and  8th  order  solutions  are  presented  and  these  results  show  this  higher 
order  boundary  conditions  are  very  promising  for  the  incompressible  flow  simulations. 
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1.  Introduction 

Today's  incompressible  flow  simulations  based  on  the  finite  difference  method  will  be  successfully 
applicable  to  the  various  practical  flows  by  using  DNS,  LES  and  RANS.  In  the  more  practical  flow 
simulations,  however,  the  enormous  computational  time  and  memory  are  necessary  in  order  to  obtain 
the  reliable  solution.  To  relax  the  restriction  on  the  limit  of  grid  resolution,  i.e.,  number  of  grid  points, 
the  higher  order  numerical  simulation  is  one  of  the  means  of  solving.  In  the  higher  order  flow 
simulation,  the  boundary  conditions,  especially  wall  boundary  conditions,  with  higher  order  of  spatial 
accuracy  are  the  more  important  problem.  One  of  boundary  treatments  is  the  use  of  one-side  finite 
difference.  In  this  case,  the  programming  is  more  complicated  because  the  various  finite  difference 
formulas  have  to  be  adopted,  and  the  compatibility  condition  for  the  pressure  usually  is  not  satisfied. 

In  this  paper,  the  new  higher  order  wall  boundary  conditions  that  introduce  the  virtual  grid  points 
inside  the  wall  according  to  the  stencil  of  higher  order  finite  difference  approximation,  are  proposed. 
In  the  present  boundary  conditions,  the  finite  difference  formulas  at  the  near  boundary  are  the  same  as 
those  at  the  inner  region.  The  square  driven  cavity  flows  are  considered  and  the  present  higher  order 
wall  boundary  conditions  are  verified. 

2.  Computational  Method 
2.1  Higher  order  flow  solver 

In  this  paper,  we  consider  the  incompressible  flows.  Then,  the  incompressible  Navier-Stokes  equations 
can  be  written  by 
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Equations  (1)  and  (2)  denote  the  continuity  equation  and  the  momentum  equations.  The  velocity 
components  and  the  pressure  are  expressed  by  ut  and  p ,  respectively.  Re  denotes  the  Reynolds  number 
defined  by  Re=UL/v ,  where  v  is  the  kinematic  viscosity.  The  equations  are  nondimensionalized  by  the 
reference  length  L ,  the  reference  velocity  U  and  the  reference  pressure  pU 2. 

The  variable  order  method  of  lines  [1,2]  is  adopted  for  solving  the  incompressible  Navier-Stokes 
equations.  In  the  method  of  lines  approach,  the  spatial  derivatives  are  discretized  by  the  appropriate 
scheme,  so  that  the  partial  differential  equations  (PDEs)  in  space  and  time  are  reduced  to  the  system  of 
ordinary  differential  equations  (ODEs)  in  time.  The  resulting  ODEs  are  integrated  by  the  Runge-Kutta 
type  time  integration  scheme. 

In  the  spatial  discretization,  the  convective  terms  are  approximated  by  the  variable  order  proper 
convective  scheme  [3],  because  of  the  consistency  of  the  discrete  continuity  equation,  the  conservation 
property,  and  the  variable  order  of  spatial  accuracy.  This  scheme  is  the  extension  of  the  proper 
convective  scheme  proposed  by  Morinishi  [4]  to  the  variable  order.  The  variable  order  proper 
convective  scheme  can  be  described  by 
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where  M  denotes  the  order  of  spatial  accuracy,  and  the  operators  in  Eq.  (3)  are  defined  by 
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where  In  this  scheme,  we  can  obtain  the  arbitrary  order  of  spatial  accuracy  by  changing 

only  one  parameter  M.  The  coefficients  c ^  and  cm>  are  the  weighting  parameters  and  Axj  denotes  the 
grid  spacing  in  the  Xj  direction. 

On  the  other  hand,  the  diffusive  terms  are  discretized  by  the  modified  differential  quadrature 
(MDQ)  method  [5].  In  the  MDQ  method,  the  spatial  derivatives  are  approximated  by  the  linear 
combination  of  the  appropriate  coefficients  and  the  function  itself.  The  second  derivative  can  be 
discretized  by 


q2u.  m/2 

y\x=  X  (X)ui\xi+m 
dxj  m=—M  /  2 


The  coefficient  ®m”(x)  is  the  second  derivatives  of  the  function  defined  by 
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The  coefficients  of  the  variable  order  proper  convective  scheme,  c ^ ,  can  be  determined  automatically 
by  using  the  MDQ  coefficients  for  the  first  derivatives.  For  example,  these  coefficients  can  be 
determined  as  follows; 
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C/  =  1  (2nd  order) 

Cj  =  9/8 ,  c3=  -1/8  (4th  order) 
cj  =  150/128,  c3=  -25/128,  c5=  3/128  (6th  order) 
d  =  1225/1024,  c3= -245/1024,  c5=  49/1024,  Ct= -5/1024  (8th  order) 

Cl  =  39690/32768,  c3=  -8820/32768,  c5=  2268/32768,  Cy= -405/32768,  c9=  35/32768  (10th 

order) 

Then  the  partial  differential  equations  in  space  and  time  are  reduced  to  the  system  of  ordinary 
differential  equations  (ODEs)  in  time.  The  resulting  ODEs  in  time  are  integrated  by  the  appropriate 
time  integration  scheme,  e.g.,  the  Runge-Kutta  type  scheme.  In  this  paper,  the  collocated  grid  system 
which  all  variables  are  defined  at  cell  center  is  adopted.  The  fractional  step  technique  [6]  is  used  for 
the  solution  procedure.  In  the  first  step,  the  fractional  step  velocity  u ,■*  is  computed  by  the  relation, 

u*=u?  +  yAtFtn  ,  (10) 


where  Ft  denotes  the  convective  and  diffusive  terms,  and  y  is  the  coefficient  determined  by  the  time 
integration  scheme.  The  superscript  n  and  *  denote  the  values  at  t  =  nAt  and  the  fractional  step  values, 
respectively.  Equation  (10)  is  solved  on  the  collocated  grid  points.  The  computed  fractional  step 

velocity  is  interpolated  into  the  staggered  locations  so  that  the  staggered  fractional  step  velocity,  uf  * , 
can  be  obtained.  By  using  Eq.  (5),  the  variable  order  interpolation  is  constructed.  The  staggered 
velocity  at  next  time  step  is  given  by  the  relation, 
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Substituting  Eq.  (11)  into  the  discrete  continuity  equation,  the  pressure  equation 


aV+7 
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can  be  obtained.  This  pressure  equation  is  solved  by  the  variable  order  SOR  method.  Finally,  the 
velocity  at  next  time  step  is  computed  by 


n+1  *  dp 

Uf  =  Uj  -  yAt 


n+1 


dxi 


(13) 


In  practice,  the  higher  order  time  integration  scheme,  e.g.,  3rd  or  4th  order  Runge-Kutta  type 
scheme,  should  be  adopted.  In  this  paper,  however,  the  forward  Euler  time  integration  scheme  (y  =  1) 
is  used,  because  we  consider  only  the  steady  flow  problem. 

2.2  Higher  order  wall  boundary  conditions 

In  the  collocated  grid  system,  the  wall  boundary  location  is  different  from  the  grid  point,  that  is,  cell 
center,  shown  in  Fig.  1.  Then,  in  order  to  estimate  the  velocities  on  the  virtual  grid  points  inside  the 
wall,  it  is  necessary  that  the  interpolated  velocities  are  equivalent  to  the  wall  (boundary)  velocities. 

Generally,  the  interpolation  is  defined  by  the  linear  combination  of  the  simple  averages  and 
coefficients. 
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where  M  and  c ^  are  the  same  ones  in  Eq.  (3)  and  t’=  21- 1 .  The  coefficients  satisfy  the  relation, 

M/2 

I ]ce>=l  (15) 

1=1 

By  using  this  relation,  the  velocity,  e.g.,  u ,  on  the  virtual  grid  point  is  specified  by 


u  -  ~(ui-r/2  +  ui+r/2)  -  uwall 


(16) 


Then,  the  interpolated  velocity  at  the  wall  boundary  satisfies. 


d  grid  system 


M/2 

uwall  =  YaCi'u 
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On  the  other  hand,  the  first  derivative  of  pressure  at  the  boundary  point  is  necessary  to  solve  the 
pressure  equation.  The  pressure  on  the  virtual  boundary  point  has  to  be  specified  by  the  first 
derivatives  of  pressure.  The  first  derivative  can  be  estimated  by  the  linear  combination  of  the  second 
order  finite  difference  with  different  grid  spacings  and  coefficients. 
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where  lv:  is  the  grid  spacing  in  the  xt  direction.  Then,  we  set  that  the  each  second  order  finite 
difference  is  equal  to  the  boundary  value, 
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so  that  the  global  Neumann  condition  can  be  satisfied, 


125 


Hidetoshi  NISHIDA 


dp 

dxt 


M/2 

Icf 

1=1 


StP_ 

Srxi 


F wall 


(20) 


where  Fwau  denotes  the  convective  and  diffusive  terms  at  the  wall  boundary. 

Moreover,  in  the  case  computing  the  second  derivative  of  velocity  at  the  boundary,  the  velocities 
on  half  grid  location  shown  in  Fig.  2  are  necessary.  In  this  paper,  these  velocities  are  determined  by 
using  the  global  conservation  [7].  The  global  conservation  is  that  the  mass  conservation  has  to  be 
specified  only  by  the  inflow  and  outflow  at  the  both  boundaries,  e.g.,  when  we  consider  the  x  direction, 
the  global  conservation  can  be  written  by 

YjAx—\i=-ui/2+uN+l/2  ,  (21) 

i=l  dx 

where  N  is  the  number  of  grid  points  in  the  x  direction.  The  subscripts  1/2  and  N+l/2  denote  the  left 
and  right  boundary  values,  respectively.  In  the  2nd  order  of  spatial  accuracy,  the  above  global 
conservation  is  satisfied  automatically,  because  of 

W&M  Boundary 


O  ;  Known  velocity  (collocated  location)  ®  :  Known  velocity  (half  grid  location) 

O  ;  Unknown  velocity  (collocated  location)  :  Unknown  velocity  (half  grid  location) 

Fig.  2  Grid  points  location 
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ox  Ax 


However,  it  is  not  satisfied  generally  in  the  higher  order  of  spatial  accuracy.  For  example,  the  global 
conservation  in  the  4th  order  is  written  by 


^  ,  du  ,4th  ,  v 
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i=l  dx 
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Equation  (23)  does  not  satisfy  Eq.  (21)  automatically.  Then,  in  order  to  coincide  Eq.  (21)  with  Eq.  (23), 
we  need  the  following  relation, 


—  (~u-l/2  ~ul/2  ~u3/2)  ~  ~ul/2  ’ 

—  (UN-1  /  2  +  UN+1  /  2  +  uN+3  /  2)  =  UN+1  /  2 
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Therefore,  the  velocity  components  at  half  grid  location  inside  the  wall,  u.1/2  and  uN+3/2 ,  can  be 
estimated.  And  in  order  to  compute  the  second  derivatives  at  the  boundary  by  usual  centered  finite 
difference,  the  velocity  components,  u.3/2  and  uN+5/2 ,  have  to  be  determined.  These  can  be  specified  by 
considering  the  continuity  equation  at  the  closest  cell  inside  the  wall  to  the  boundary.  In  the  case  of 
Fig.  2,  the  following  relation  can  be  obtained  from  the  continuity  equation  on  ( 0,j ). 

du  .  C;  ,  c3  ,  v  5v  , 

~dx  °J  ~  ~Ax  (~U~1/2’j  +Ul/2’j)  +  YAx  ~U~3/2j  +U3/2’j)  ~~~dy  0,J  25 

The  right  hand  side  of  Eq.  (25)  is  known  by  using  the  interpolated  values.  Then,  the  velocity,  u.3/2  ,  can 
be  estimated.  Similarly,  the  higher  order  wall  boundary  conditions  than  the  4th  order  of  spatial 
accuracy  can  be  formulated. 

3.  Computational  Results 

In  this  paper,  the  square  driven  cavity  flow  problem  is  considered.  The  computational  conditions  are 
that  the  number  of  grid  points  is  41x41  (uniform  grid),  Reynolds  number  Re=1000  and  5000 ,  and  the 
convergence  criteria  is  L2-residual<10'6  for  the  Navier-Stokes  and  pressure  equations. 


(a)  2nd  order  solution 


(b)  4th  order  solution  (c)  6th  order  solution  (d)  8th  order  solution 
Fig.  3.  Pressure  coefficient  distributions  ( Re=1000 ) 


Velocity  (u) 

Fig.  4.  Velocity  profiles  along  the  center  lines 
( Re=1000 ,  41x41  grid  points) 
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Fig.  5.  Velocity  profiles  near  the  boundaries 
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Figures  3  and  4  show  the  pressure  coefficient  distributions  and  the  velocity  profiles  along  the  center 
lines  with  the  Reynolds  number  Re=1000  and,  the  2nd,  4th,  6th,  and  8th  order  of  spatial  accuracy, 
respectively.  The  pressure  coefficient  is  defined  by 


Cp  = 


P~Po 
1  /  2pU2 


(26) 


where  p0  is  the  pressure  at  center  of  bottom  cavity  wall  and  U  denotes  the  moving  wall  velocity.  In  Fig. 
4,  the  cross  symbol  denotes  the  result  with  129x129  grid  points  obtained  by  Ghia  et  al.  [8].  The 
velocity  profiles  near  the  boundaries  are  plotted  in  Fig.  5.  The  circular,  triangular,  inverse  triangular, 
and  square  symbols  denote  the  2nd,  4th,  6th,  and  8th  order  solutions,  respectively.  In  this  Reynolds 
number,  the  higher  order  solutions  than  the  second  order  are  almost  the  same.  In  comparison  with  the 
Ghia's  solution,  these  higher  order  velocity  profiles  are  clearly  improved. 


Figures  6,  7,  and  8  show  the  results  in  Reynolds  number  Re=5000 ,  that  is,  the  pressure  coefficient 
distributions,  velocity  profiles  along  the  center  lines,  and  velocity  profiles  near  the  boundaries.  In 
contrast  with  the  previous  case,  the  independent  solution  of  spatial  accuracy  can  not  be  obtained. 


(a)  2nd  order  solution  (b)  4th  order  solution  (c)  6th  order  solution  (d)  8th  order  solution 

Fig.  6.  Pressure  coefficient  distributions  ( Re=5000 ) 


tt  i  i  Velocity  (u) 

However,  as  the  spatial  accuracy  becomes  higher,  Fig  8  Veiocity  proflies  near  the  boundaries 

the  primary  vortex  becomes  stronger  and  the 


128 


Higher  Order  Wall  Boundary  Conditions 


velocity  profiles  become  closer  to  the  Ghia's  solution.  Especially,  the  velocity  profiles  near  the 
boundaries  are  different  between  the  orders  of  spatial  accuracy.  These  distributions  show  the  tendency 
to  converge  at  the  highest  order  solution.  Then,  the  present  higher  order  wall  boundary  conditions  are 
effective  to  improve  the  solutions. 

4.  Concluding  Remarks 

In  this  paper,  the  higher  order  wall  boundary  conditions  for  incompressible  flow  simulation  are 
presented  on  the  collocated  grid  system.  The  validation  is  performed  in  the  steady  square  driven  cavity 
flow  problems,  i.e.,  Re=1000  and  5000  with  41x41  grid  points,  by  using  the  variable  order  method  of 
lines. 

In  the  case  of  the  Reynolds  number  Re=1000 ,  the  independent  solution  of  spatial  accuracy  can  be 
obtained.  On  the  other  hand,  as  the  order  of  spatial  accuracy  becomes  higher,  the  primary  vortex 
becomes  stronger  and  the  velocity  components  near  the  boundaries  are  clearly  improved  in  the  case  of 
the  Reynolds  number  Re=5000. 

Then,  it  is  concluded  that  the  present  higher  order  wall  boundary  conditions  are  more  available  for 
the  incompressible  flow  simulations.  In  this  paper,  the  higher  order  boundary  conditions  are 
constructed  on  the  collocated  grid  system,  but  it  is  possible  to  build  the  almost  the  same  boundary 
conditions  on  the  staggered  grid  system.  Therefore,  this  idea  of  higher  order  boundary  conditions  is 
very  promising  for  the  higher  order  incompressible  flow  simulations. 
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